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We consider steady gravity-driven flow of a thin layer of viscous fluid over a curved 
substrate. The substrate has topographical variations ('bumps') on a large scale com- 
pared to the layer thickness. Using lubrication theory, we find the velocity field in 
generalized curvilinear coordinates. We correct the velocity field so as to satisfy 
kinematic constraints, which is essential to avoid particles escaping the fluid when 
computing their trajectories. We then investigate the particle transport properties 
of flows over substrates with translational symmetry, where chaotic motion is pre- 
cluded. The existence of trapped and untrapped trajectories leads to complicated 
transport properties even for this simple case. For more general substrate shapes, 
the trajectories chaotically jump between trapped and untrapped motions. 
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I. INTRODUCTION 

It is often the case in applications that viscous fluid flows in a thin layer at low Reynolds 
number, which allows lubrication (or long-scale, also long- wave) theory to be used. 1 The 
resulting decoupling between the long-scale and the 'thin' direction greatly simplifies the 
mathematical solution of the flow. There is essentially only one equation to be solved, 
for the thickness of the fluid as a function of position. This equation arises from mass 
conservation appropriately averaged over the thin coordinate. The velocity is then deduced 
from the thickness. 

The mass-conservation equation to be solved involves only the long-scale coordinates. The 
third coordinate is 'slaved' to the horizontal ones by the effect of viscosity. In that sense one 
usually thinks of thin-layer flows as 'two-dimensional,' in the sense that vertical variations 
(i.e., in the thin direction) exist but are subsumed in the theory and that the vertical velocity 
field is small. However, if one is interested in the transport properties of such flows, any 
amount of vertical variation is crucial. We will show here that even small variations can lead 
to chaotic behavior of particle trajectories, known as chaotic advection, 2-4 greatly enhancing 
the transport properties of the flow. Indeed, if the flow were truly two-dimensional and 
steady, chaotic transport would be impossible since two-dimensional autonomous dynamical 
systems cannot exhibit chaos. 5 

As far as we know, the transport properties of thin-layer flows have not been studied, 
beyond plotting streamlines in two-dimensional cases, probably because they were assumed 
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too simple to be important. But understanding these properties could be important if, 
for instance, a compositionally-uniform coating is required in an industrial application. A 
related area that has received much attention lately is chaotic advection in microchannels. 6 
There periodic grooved patterns at the bottom of the channel cause swirl in the flow that 
leads to chaos. Here we study a different regime where the bottom substrate has larger 
topographical variations, the is fluid thinner, and the sidewalls less important. 

We will show here that the shape of the curved substrate is crucial. Though we investigate 
this for the case of steady gravity-driven flow with a free surface, the conclusions reached 
should qualitatively apply to different thin-flow situations such as those driven by surface 
tension gradients. Constant surface tension is included in the derivation, but left out of 
most of the examples since it does not affect the qualitative results. 

There are several derivations of models of thin fluid flow over curved substrates in the 
literature. 7-12 Our approach follows that of Roy et al. 8 (hereafter RRS), in that the vertical 
(thin-direction) coordinate is orthogonal to the long-scale coordinates, but differs in that 
the long-scale coordinates are not necessarily orthogonal to each other. This is because 
even though orthogonal coordinates can always be found, 13 doing so for a general substrate 
(other than planes, cylinders, and spheres) requires solving differential equations numerically. 
Working in generalized coordinates is not much more complicated if differential-geometric 
notation is used, 14-16 and allows a simple parametrization of the substrate in terms of a 
height function. 

Another motivation for re-deriving the model of RRS 8 is that we need to include correction 
terms to the mass-conservation equation and vertical velocity to ensure that the kinematic 
boundary condition at the free surface is satisfied exactly, up to numerical accuracy. In the 
RRS model, the kinematic boundary condition is satisfied to second order in the thickness. 
This is fine for most applications, but for particle advection it is catastrophic: in numerical 
simulations, individual particles will gleefully escape the top surface if they come anywhere 
near it, since even a second-order error is enough to give a small component of the velocity 
field normal to the surface. Rather than remedy this in some ad-hoc manner (by reflecting 
the particles back in, for instance), we instead include some judiciously-chosen second-order 
terms that guarantee that the kinematic boundary condition at the top surface is satisfied. 

The outline of the paper is as follows. Sections ILTV comprise the derivation of the model 
equations. At each stage, we carefully record where the thinness of the layer comes in as 
an approximation. Since we will be working with substrates of complicated shapes, the first 
part of of our paper (Section II) is devoted to establishing a convenient coordinate system 
for our problem. Then we tackle mass conservation in those coordinates in Section III. 
Section IV introduces the dynamical equations of motion for Stokes flow. In that section 
we also rescale the variables (Section IV A) and solve for the velocity field (Section IV B) to 
first order in the thickness. This gives the mass flux vector, whose divergence must vanish 
for steady flow. We include second-order terms to the mass flux to impose the kinematic 
boundary condition at the free surface, which is essential for particle advection as discussed 
above. 

In Section V we discuss solutions of the governing PDE for the fluid thickness. The 
simplest type of solution involves substrates with a translational symmetry, for which the 
steady flow can be found analytically at leading order. However, solving for the flow at first 
order, or for substrates without such a symmetry, we must resort to numerical solutions. 

Section VI is concerned with fluid particle advection and is the central results section 
of the paper. We examine the variety of particle trajectories by making Poincare surfaces 



3 



of section, for increasingly complicated substrates. As expected, the transport properties 
also become more complicated, exhibiting a range of behavior from trapped to chaotic 
trajectories. We offer some concluding remarks in Section VII. 



II. COORDINATE SYSTEM 

A. Separating the Thin Direction 

In our problem, fluid motion occurs over a curved substrate of arbitrary shape. The 
direction normal to the substrate is special in that it defines the direction in which the fluid 
layer is assumed 'thin.' Hence, it is convenient to locate a point r in the fluid as 

r(x\ x\ y) = X{x\x 2 ) + y e 3 (x\ x 2 ) (II.l) 

where X(x 1 ,x 2 ) is the location of the substrate, e 3 is a unit vector normal to the substrate, 
and y is the perpendicular distance from r to the substrate. The coordinates x 1 and x 2 are 
substrate coordinates used to localize points on the substrate. For example, in Section II B 
we will use the Monge parametrization, X = [x 1 x 2 /(x 1 ,^ 2 )) , where / gives the height 
of the substrate. 

The tangent vectors to the substrate are 

e a = d a X (II.2) 

where d a := d/dx a . The coordinate vectors associated with the coordinate system are found 
from (II.l), 

e a = d a r = e a -y Kj ep , — = e 3 , (II.3) 

where we assume the summation of repeated indices, and Greek indices only take the value 1 
or 2. We have also defined the curvature tensor 1K Q /3 from 

d a e 3 = -Kje , (II.4) 

since changes in the unit vector e 3 must be perpendicular to e 3 , and e 3 • e a = 0. Note that 
the e a are not necessarily orthogonal or normalised. 

We adopt the convention that quantities with a tilde are evaluated in the 'bulk' (away 
from the substrate), and thus depend on y, whilst those without the tilde are 'substrate' 
quantities and do not depend on y. Thus, e^x 1 , x 2 , 0) = e Q ,(a; 1 , x 2 ). Starting in Section III, 
we will also denote quantities evaluated on the free surface at y = i] by an overbar. 

The three-dimensional metric tensor Q ab has components 



where 



da/3 = e a ■ ep = G a/3 , g a3 = e a ■ e 3 = , g 33 = e 3 • e 3 = 1 , 



G af3 := e a - ep = G a p - 2y K af3 + {y} 2 K a 7 K 7/3 , 
^ap '■= e a - ep , 



(II.5) 



and we have written the symmetric tensor K a p = K a 7 Gp^ in the usual manner of lowering in- 
dices with the metric tensor. 14-16 The three-dimensional metric tensor is thus block-diagonal, 
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and the y coordinate is unstretched compared to the Cartesian coordinate system. It mea- 
sures the true perpendicular distance from the substrate to a point in the fluid. Note that it 
is important that indices be raised or lowered with the appropriate metric: G a/ g for substrate 
quantities, Q a b or G a/ g for bulk quantities. The tensor G a p(x 1 , x 2 , y) can be used to raise or 
lower indices, even though it depends on three coordinates but only has two-dimensional in- 
dices, because of the block-diagonality of the metric q^. Note also that in Eq. (II. 5) we have 
written {y} 2 rather than y 2 , since the second form could be misinterpreted as the second 
contravariant component of y (which is nonsense in this case). We shall follow this practice 
throughout the paper, but note that negative powers are written as y~ 2 , since then there is 
no confusion. 

Given the substrate vectors e a , it is easy to solve for the covectors e a , which are such 
that e a ■ ep = 5p a . Then the bulk covectors are 

e a = e a + y^ a e (3 + 0{{y} 2 ), (II.6) 

to first order in y, which can be checked by taking their dot product with the coordinate 
vectors (II. 3) to see that e a ■ ep = 8p a + 0({y} 2 ). We will not need the covectors to higher 
order in y. From (II. 6), we find the inverse metrics, 



(II.7) 



QaP = ga.g/3 = G «P + 2 y + ({y} 2 ) , 

G a/3 = e a -e /3 , 

where again the higher-order terms will not be needed. 
A crucial quantity is the determinant of the metric, 

w = (det G afS ) 1/2 = Cuw, w= (det G a/3 ) 1/2 , (II.8) 

where 

u = l-Ky + ${y} 2 , (II.9) 

and 

k = K a a mean curvature (II. 10a) 

9 = det Kf/ Gaussian curvature. (II. 10b) 

The area element at fixed y is given by u w dx 1 dx 2 , and the volume element by Cj w dx 1 dx 2 dy, 
since the y coordinate is unstretched with respect to the Cartesian system. The volume per 
unit substrate area of a fluid layer of thickness rj is thus 

c-£5^-,-W,}-W. (n.n, 

The quantity u also tells us when the coordinate system ceases to be valid: if Cj — 0, then 
the area element becomes singular. Hence, the valid range for y is 

< y < {max^x, k 2 , 0+)}" 1 , fc 1>2 := \(k ± ^ {k} 2 - 4S) 

where k\ and k 2 are the principal curvatures. 13 If both principal curvatures are negative then 
the positive range in y is unrestricted (convex region); otherwise y is limited by the largest 
positive principal curvature. 

We now have all the geometrical equipment we need to solve fluid equations on the curved 
substrate. Note that all the results of this section are exact, except for the covectors (II. 6) 
and inverse metric (II. 7), which are valid to first order in y. 
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B. Substrate Coordinates 



For most applications in the literature of thin films, orthonormal coordinates have been 
the coordinates of choice. This is because the main substrate shapes that have been treated 
are planes, cylinders, and spheres, where orthonormal coordinates are readily available. For 
a general substrate shape, orthonormal coordinates are difficult to construct and require nu- 
merical integration. Singularities (umbilics) also cause problems. 13 For our application — flow 
down a curved substrate — the Monge representation of a surface 17 is the most convenient. 

The Monge representation is a glorified name for a parametrization of the substrate by 

X(x\x 2 ) = (x 1 x 2 f{x\x 2 )f (11.12) 

in three-dimensional Cartesian space. Following standard notation, 17 we define 

p = dj, q = d 2 f, (II. 13a) 

r = d 1 d 1 f, s = d 1 d 2 f, t = d 2 d 2 f. (II.13b) 

The unnormalized, nonorthogonal tangents e\ and e 2 are 

e 1 = d 1 X= (1 pf, e 2 = d 2 X=(0 1 qf , 

and their normalised cross product gives the normal to the substrate, 

e 3 = - {-p -g 1) . 

The corresponding covectors are 

e 1 = j\~ 2 ((1 + {q} 2 ) -pq pf , e 2 = * (-pq (1 + {p} 2 ) qf 
and e 3 is its own covector. The metric tensor of the substrate and its inverse are 

(11.14) 

with determinant 

,V2 



«; = (det G Q/3 ) = v/T+MW- 

The curvature tensor is 

/?i . 1 ^(1 + {g} 2 )r — pgs (1 + {p} 2 )s — pgr 



{Ka }: ~ M 3 V(i + W 2 )s-^ (i + M 2 )t-wsJ (IL15) 



which gives mean and Gaussian curvatures 



K = K a a = j^ ((1 + {q} 2 )r - 2pqs + (1 + {p} 2 )t) , 
5 = detK J = -±-(rt-{s} 2 ). 
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We write the normalized gravity vector as g = (sin 9 cos sin 9 sin — cos9) T , so that 
the inclination angle 9 is zero for gravity pointing downwards, and for e (— 7r/2,7r/2) 
positive # induces flow in the positive x 1 direction. Then we have the components 

g] = g ■ e 1 = — (p cos 9 + pq sin 9 sin — (1 + {g} 2 ) sin 9 cos 0) / {w} 2 , 

g 2 = g ■ e 2 = — (q cos 9 + pq sin cos — (1 + {p} 2 ) sin sin 0) / {w} 2 , (11.16) 

g y = g ■ e 3 = — (cos + p sin cos + g sin sin 0) /iy . 

The specific parametrization of the substrate introduced in this section will not be needed 
in the derivation of the equations of motion (Sections III and IV), only in their solution 
(Section VB). Hence, a different parametrization could be used if called for by the geometry 
of the substrate. For instance, flow down a curved filament is better parametrized by 
cylindrical coordinates, or if the substrate has overhangs (making / multivalued) coordinates 
based on arc length are better. 



III. MASS CONSERVATION 

We remind the reader of our notational convention: 

1. Quantities with a tilde (e.g., e a and w) are evaluated between the substrate and the 
free surface (in the 'bulk') and are functions of (x 1 , x 2 , y, r). 

2. Quantities with an overbar (e.g., e a and w) are evaluated on the free sur- 
face y = ^(x 1 , x 2 , r) and are functions of (x l , x 2 , r). 

3. 'Bare-headed' quantities (e.g., e a and w) are evaluated on the substrate y = and are 
functions of (a; 1 , a; 2 ), or they are quantities that do not depend on y at all (e.g., e 3 ). 

We introduce a velocity field 

u = u a e a + v e 3 , (HI- 1) 

associated with some time-dependent incompressible flow. Here we depart from RRS 8 in two 
ways: first, the e a are not unit vectors, so the velocity components u a are scaled differently; 
second, we use e a and not e a , so that the velocity field is expressed in terms of the bulk 
coordinate vectors rather than the substrate tangent vectors. This is slightly more natural 
in the present context, but is not crucial. 

Mass conservation is imposed via the divergence-free condition, V ■ u = 0; in terms of 
the coordinates in section II, this is 

d 

d a (wu a ) + — (wv) = 0. (III.2) 

We integrate this from to r\ and use the no-throughflow condition v(x l , x 2 , 0) = to get 

w v — — [ V d a (wu a ) dy (III.3) 
Jo 

= wu a d aV -d a [ (wu a )dy. (III.4) 
Jo 
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Now we use the kinematic boundary condition at the top surface, 

d T r] + u a d a r/ = v, (III. 5) 

to find 

rv 

wd T i] = —d a / (wu a ) dy. (III. 6) 

Jo 

We define the mass flux vector 

q a {x\x\y) := f 'cou a dy, <F(x\x 2 ) = P(x\i?,ri) (III.7) 
Jo 

which allows us to rewrite (III. 6) as 

w d T 7] = -d a (w q a ) . 

We divide through by w, 

ud T r] = -V a q a (III.8) 

where the covariant divergence V a q a := iu -1 £? a (it; q a ), with V a denoting the covariant deriva- 
tive with respect to x a using the substrate metric G Q/3 . 14-16 In Eq. (III. 8), Co defined by 
Eq. (II. 9) is evaluated at y — rj, and r only enters Co through 77. Hence, we can rewrite (III. 8) 

as 

^ + W a q a = 0, (III.9) 

where ( is the volume per unit substrate area given by Eq. (11.11). Equation (III. 9) thus 
embodies mass conservation. 

Note that there were no assumptions on the thinness of the layer in this section. 



IV. DYNAMICAL EQUATIONS OF MOTION 



Sections III involved only kinematic considerations, namely mass conservation and the 
kinematic boundary conditions at the free surface and bottom substrate. Hence, everything 
we have done so far applies to any incompressible flow. Now we will restrict to the problem 
of a viscous (Stokes) flow driven by gravity, which satisfies the equation 



Au = Vp - g, 



(IV.l) 



where p is the pressure and g is a unit vector in the direction of gravity. The velocity satisfies 
the boundary conditions 



u 

t a • T • h 
—p + h ■ T • n 



d T rj + u a d a r] 



where 



at y = no-slip at substrate (IV. 2a) 

at y = i] tangential stresses at free surface (IV. 2b) 

(tk sw { at y = i] normal stress at free surface (IV. 2c) 

v at y = i] kinematic condition at free surface (IV. 2d) 

T := Vu + (V«) T 
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is the deviatoric stress, n is the unit normal to the surface, t a are tangents to the surface, 
and K snr f is the mean curvature of the surface. The stress conditions assumes that the gas 
above the surface has negligible viscosity and density. 

We have absorbed the constant density, viscosity, and strength of gravity into the ve- 
locity and pressure, so u stands for (uu/g), p stands for (p/gp), and the surface tension a 
(introduced below) stands for (a/gp). This means that the velocity in (IV. 1) has units of 
squared length, the pressure has units of length, and the surface tension has units of squared 
length. 

In the coordinates of Section II, the location of the surface is 

f(x 1 , x 2 , t) = X(x 1 , x 2 ) + rj(x l ,x 2 , r) e^x 1 , x 2 ) 

and the surface tangents i a are then given by 

i a = d a f = e a + d a r] e 3 . (IV. 3) 

The unit normal h is obtained by taking the cross product of the tangents and normalizing. 

In the coordinates introduced in Section II, with the velocity field (III. 1) , The Stokes 
Eq. (IV. 1) becomes 

i jd Q (w + |- (w | («t e 7 + v e 3 ) = e a d a p + e 3 ^ - gf e« - & e 3 , (IV A) 

where we wrote 

9 = 9s e * + 9y es 

with components given by (11.16). Note that, unlike the velocity field (III. 1) , we have not 
written the gravity vector in terms of the bulk (tilde) coordinate vectors. This is to avoid y 
dependence in g£ and g y . 



A. Small-parameter Rescaling 

Following RRS, 8 we rescale the variables to account for the thinness of the fluid layer. 
We do this not by rescaling y, but by rescaling x a : we treat the variations in the substrate 
as being large-scale, whilst the depth of the layer is of order unity. Hence, we introduce the 
rescalings 

x — e x , v — e v , p = e p , a — e o ■ , 

where e is a small parameter; all other variables are left unsealed. Following standard 
practice, we immediately drop the * superscripts and deal only with the rescaled variables. 
The amplitude of the variations of the substrate are taken to be such that the vectors e a and 
hence the substrate metric G are of order unity. All terms involving the curvature tensor IK 
are of order e, because of the d/dx a in the definition (II. 4). 

The continuity equation (III. 2) is unchanged by the rescalings, and neither is the kine- 
matic boundary condition (III. 5). The dynamical equation (IV. 4) becomes 

| { {e} 2 d a (w G^fy) + |- (w | {v? e 7 + e v e 3 ) = e a d a p + e' 1 ^ ^ - £ e a - g y e 3 , 

(IV.5) 
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where w = w(l —e ny) + 0({e} 2 ), e a = e a — ey K a @ ep, and e a = e a +eyK a/3 e /3 + 0({e} 2 ). 
As for the boundary conditions (IV. 2), they simplify to 

u a = d = at y = ; (IV.6a) 

du a 

__ =0 + O({e} 2 ) aty = V ; (IV.6b) 

p = -a K surf + ({e} 2 ) at y = rj , (IV.6c) 

where 

«Wf = « + £ V K 2 + zAr] + ({e} 2 ) 
with Ar) = V a V a r] = w~ 1 d a (w G a ^dprj) the covariant Laplacian of rj, and 

We shall not need the higher-order terms in (IV.6). 

B. Velocity Field 

We now proceed in the usual manner and expand the fields as power series in e: 

= «(0) + e ufi) + • • • , 
P = P(o) +£p(i) + •••, 

where we have left out v since it is not needed in our development until later (and it will 
not appear as an asymptotic series). At order e°, the velocity field and pressure satisfy 



d 2 u a 



(o) _ na~ ~ a d P(o) 



d{y} 2 *™ " ' dy 
where <9 a = G a P dp. These are readily integrated to give 

«fo) = ^fo) " 2/7), p ( o) = -<r«, (IV.7) 

where 

Afo = -§(£ + a 0"k) (IV.8) 

and the boundary conditions (IV.6) have been applied. 
At order e, 

d 2 u a dvP F)~ 

0£yf 2 = («V + 2 V) - + + 2y K^^ ( o) - y yf K," , -|H = ^ , 

with solution 

fifi) = Af 1} y (y - 21,) + B ( V ({y} 2 - 3{ij} 2 ) (IV.9) 

where 

= \ {9>V + 3yf - 9yd a V + a {v (k0°k + 2 K/c^) - d a (K 2 r] + Ar))}) , 

(IV. 10a) 

= "I + <r^«) (k V + 4 V) • ( IV - 10b ) 
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To simplify these expressions, we made use of 

d a g v = d a (g- e 3 ) = -g • K/e^ = -g% K afS , 

since g is a constant vector. 

Now that we have u a to first order in e, we can insert it in the definition of the mass flux 
vector (III.7), 

q a = [ ! uu a dy= f\l-eKy)(uf 0) +euf 1) )dy + O({e} 2 ), (IV.ll) 
Jo Jo 

which gives 

r = Cv + Crf (IV.12) 

with 

= |M 3 {g: - e (« V* + \ V) v + eg y d a v } 

+ {4 2 iio M 4 « flf ( 9k V + 11 V) ~ 25 9v d a v} + ({e} 2 ) , (IV.13) 



+ i £ } 2 lk ° ivV K i 9r l K afXK - ^r]K p a d p K - 25 d a (K 2 r] + Ar])} + ({e} 2 ) , (IV.14) 

where we have kept some second-order terms but neglected others: this will be explained 
when we find v below and is related to the exact preservation of the kinematic constraints of 
the problem. Equation (III. 9) with the flux (IV.12) is the dynamical equation that governs 
the evolution of rj. To order e, it is the same equation as that found by RRS, 8 but in 
nonorthogonal coordinates rather than orthogonal coordinates. This is a straightforward 
generalization, but it is illuminating to rederive the equations in this setting, and allows a 
thorough introduction to all the notation. We will also need the quadratic terms in (IV.12) 
below, and these are not present in RRS's formulation: this is not a failure of RRS but 
rather is required by the specific application we have in mind, namely particle advection. 

In order to do particle advection in Section VI, we will need the vertical velocity field, v. 
It is obtained by integrating the continuity equation (III. 2), 

v = -L J* d a (wu a )dy (IV.15) 

where the integral's lower bound enforces the no throughflow boundary condition at y — 0. 
The kinematic boundary condition (III. 5) is satisfied by (IV.15) at y = rj, since 

i r i r i 

v = — - d a (wu a )dy = — -d a I (w u a ) dy + — (w u a )d a r) 
w Jo w Jo w 

= -^d a (wq a )+u a d a r]. (IV. 16) 

The first term on the right-hand side of (IV. 16) vanishes because of (III. 9), and the second 
gives the kinematic boundary condition (III. 5). But observe that to recover this boundary 
condition we need factors of id to cancel; hence, a factor of w must appear multiplicatively 
in the integrand if we are to satisfy the kinematic boundary condition exactly (i.e., to all 
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orders in e). This is the reason for keeping some second-order terms in (IV. 11) and (IV. 12). 
Note that this way the divergence-free condition (III. 2) is also satisfied exactly. 

The easiest way to get the y coefficients for v is to combine (IV. 7) and (IV. 9) into 

u" = fi « 0) + eVfa = A a y (y - 2rj) + B a y ({y} 2 - 3{ V } 2 ) , (IV.17) 

valid to first order in e, where A a = A^L+eA?^ and B a = eB?^. Then, using w = w(l—EKy), 
we find from (IV. 15) 

v = -r^— (a + by + c {y} 2 + d {y} 3 ) (IV.18) 

where 

a = ±V a (2rjA a + 3{ V } 2 B a ) , c = |V Q {enA a - B a ) , 

b = -|V a ((1 + 2eKr])A a + 3eK{r]} 2 B a ) , d = |V Q (ekB*) , 

and the covariant divergence is defined after Eq. (III. 8). The expression (IV.18) for v is 
exact, given the first-order expressions w = w(l — eny) and (IV.17) for u a . 



V. SOLUTION OF THE EQUATIONS 

We now seek steady solutions of Eq. (III. 9) with the mass flux vector (IV. 12). The 
unknown function is the thickness field 77 (x 1 , x 2 ). We characterize the shape of the substrate 
using the parametrization in Section II B, and assume the substrate is periodic in x 1 and x 2 
with unit period. In dimensionless form, the independent input parameters of the problem 
are thus 

1. The inclination angle 9 and rotation angle of the gravity vector [Eq. (11.16)]; 

2. The thickness of the film at the corner, 77(0, 0) = e; 

3. The shape of the substrate, including its height, specified with /(a; 1 ,^ 2 ); 

4. The surface tension a. 

The parameters 1 determines the speed and direction of the flow: the fluid is static for 9 = 0. 
Parameter 2, the thickness of the film at the corner, effectively determines e. We shall thus 
drop e from the equations for the remainder of the paper, instead introducing the small 
parameter through i](0, 0) = e. Parameter 3 is infinitely rich: even for a given substrate 
shape, we can still vary its height by scaling / by a constant. We shall limit ourselves to 
substrates of relatively simple harmonic shapes. 



A. Substrates with Translational Symmetry 

Consider first the case where the substrate has a translational symmetry. We may then 
orient our coordinates such that / is only a function of x 2 . From Eq. (11.13), we then 
have p = s = t = 0. Since nothing depends on x 1 , we seek solutions of the form i](x 2 ). 
Equation (III. 9) then reduces to 

d 2 (wq 2 ) = =^ wq 2 = const. (V.l) 
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FIG. 1: Surface thickness tj(x 2 ) for a substrate with translational symmetry f(x 2 ) = 0.1 sin(2-7rx 2 ), 
for e = 0.02, = 0.7, (p = 1.2, and a = 0. The solid line is the zeroth-order analytic solution (V.2), 
and the dashed line is the first-order numerical solution. 

since d^wq 1 ) vanishes. If we keep only the e° terms in (IV. 12), we can solve (V.l) analyti- 
cally, 

r](x 2 ) = {cw{gl + aw' 2 d 2 K)Y l/Z (V.2) 

where the constant c is adjusted to satisfy the boundary condition 7/(0) = e, and we have 
used G 22 = w~ 2 . The quantity in parentheses in (V.2) can become singular or negative, an 
indication that the higher-order terms can no longer be neglected. In that case ^(wg 2 ) = 
still holds but we then have to solve an ODE for r](x 2 ). We shall not do this here and instead 
solve the general PDE in the next section. For comparison, the zeroth-order solution (V.2) 
is plotted along with the numerical first-order solution for a sinusoidal substrate and corner 
thickness s = 0.02. The two agree well since the thickness is small. 

Note that even though the velocity field depends only on two variables (x 2 and y), the 
flow is actually three-dimensional, since A 1 and B l in (IV. 17) do not in general vanish, 
except in the special case = tt/2. We shall see in Section VI that even this simple solution 
supports rather complicated trajectories. 

B. General Numerical Solution 

In general, Eq. (III. 9) with the mass flux vector (IV. 12) must be solved numerically. We 
use a fourth-order finite-difference scheme together with Newton-Kantorovich iteration 18 
to converge to the steady solution. Since the N-K scheme employs a linearization of the 
nonlinear equation to seek an improved guess, it allows us to check the linear stability of the 
flow for free, and all the flows presented here are stable. In Figure 2 we show the thickness 
field of a typical solution, for the substrate 

f{x\ x 2 ) = 0.05 (sin(27rx 1 ) sin(2vrx 2 ) + 0.2 sin(4vrx 2 )) (V.3) 
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(a) (b) 



FIG. 2: The thickness field r](x l , x 2 ) for the substrate shape given by Eq. (V.3), with e = 0.06, 9 = 
0.1, <j) = 0, and a = 0. (a) Contour plot: the flow is towards the right; (b) Surface plot: the flow 
is towards the upper-right. The vertical scale is exaggerated. 

with e = 0.06, 9 — 0.1, <fi — 0, and a = 0. As mentioned in the introduction, we will set 
the surface tension to zero. The solutions presented in this paper typically have resolutions 
of 100 x 100 gridpoints, which is well-resolved enough. The acid test to see if we have enough 
numerical resolution is whether particle orbits 'spiral in' to attractors, which is precluded 
if the flow is incompressible. For example, the orbits in the regular islands of Fig. 5(a) are 
nice and closed over many periods. 

When doing particle advection in Section VI the thickness field is first resampled on 
a grid of 1024 x 1024 using Fourier interpolation, and then the velocity field is evaluated 
off-gridpoints using linear interpolation. This has the advantage that the hardest part (the 
resampling) only needs to be done once at the beginning and stored. This is cheaper, for 
instance, than using a higher-order interpolation scheme on the coarser grid. 

VI. FLUID PARTICLE TRAJECTORIES 

To investigate the transport properties of the flow, we solve the advection equation 

d T r(T) = u(r(r)), r(0) = r , (VI. 1) 

for a steady velocity field u of the form (III. 1) obtained numerically as described in Sec- 
tion V. Instead of converting our velocity field back to Cartesian coordinates, we use the 
definition (II. 1) of r(r) and take its time derivative, assuming that the coordinates are a 
function of time, and use the chain rule: 

d T r = d a r d T x a + d y r d T y = d T x a e a + d T y e 3 , 



so that (VI. 1) becomes 



d T x a = u a , 



d T y = v 



(VI.2) 
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(a) (b) 

FIG. 3: For the same parameters as in Fig. 1: (a) Poincare section at x 1 = 0; (b) Four typical 
trajectories at different depths y. The trajectories have the same overall slope. 

We can thus follow particle trajectories in out curved coordinates using (VI. 2). We will 
now do so for different types of substrates: those possessing translational symmetry (Sec- 
tion VIA) and those leading to chaotic trajectories (Section VI B). In all cases, the parti- 
cle trajectories are computed using an adaptive embedded fourth/fifth-order Runge-Kutta 
scheme 19,20 and the velocity field is interpolated as described in Section VB. 

A. Substrates with Translational Symmetry 

For substrates with translational symmetry, as described in Section V A, the system (VI. 2) 

is 

d T x 1 = u 1 (x 2 ,y), d T x 2 = u 2 {x 2 ,y), d T y = v{x 2 ,y). (VI.3) 

We may immediately write x 1 = tt 1 (a; 2 , y) r + xj, so that (VI.3) is effectively an autonomous 
two-dimensional system, which is solvable by quadratures and cannot exhibit chaos. 5 The 
system admits a streamfunction ip(x 2 ,y) such that 

u 2 = w -1 dytp, v = —w^ 1 d2ip 

from which it is easily seen that ip = wq 2 , defined in (III. 7). If the zeroth-order terms 
dominate everywhere, the resulting velocity profile is quadratic in y ('half-Poiseuille') and the 
flow in the x 2 -y plane cannot have a stagnation point inside the fluid, except at y = 0. Thus, 
a Poincare section of trajectories, made by recording intersections with the plane x 1 = 0, 
must look like Fig. 3(a), where every trajectory crosses the domain in x 2 by going over 
furrows and crests. In Fig. 3(b) we see that four typical trajectories have roughly the same 
overall 'slope.' This leads to poor transport properties, since particles initially near each 
other remain near. This is due to the thinness of the layer which causes the zeroth-order 
velocity to dominate. If we use the zeroth order horizontal velocity (IV. 7), we have 

^(0)/^(0) = ^(o)M(o) 
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FIG. 4: For the same parameters as in Fig. 1 but with 9 = 0.6, (f> = 0.8: (a) Poincare section 
at x 1 = 0, showing no apparent chaos; (b) The types of trajectories seen in (a): the bottom 
trajectory is trapped in a furrow, whilst the top two trajectories moves rapidly over crests, and more 
slowly over furrows. The three trajectories start at different depths y. Four typical trajectories, 
showing a trapped trajectory and three trajectories that undergo regular 'near-trapping' events. 
Note that the overall slopes of the trajectories are different, and are not monotonic with depth. 

independent of the depth y. This implies that when the layer is very thin the trajectories 
in Fig. 3(b) are independent of depth. (This is true for any substrate, not just those with 
translational symmetry.) First-order corrections terms are required in order to have more 
complicated dynamics. 

If the inclination angle 9 of the substrate is decreased, or <p is changed so that the furrows 
present a more acute angle to the incoming flow, then a situation as in Fig. 4 develops. Here 
some of the fluid cannot make it over the crests and is trapped in a recirculating region. 
The trajectories that are not trapped go very quickly when over crests, as indicated by the 
paucity of points in the shallow regions of the Poincare section Fig. 4(a), and linger over 
furrows, where the fluid is deeper. This is of course a consequence of mass conservation, 
as the fluid is forced to flow more rapidly in the shallows. Three representative trajectories 
are depicted in Figure 4(b), where the change in speed in the untrapped trajectory (top) is 
reflected in the slope of the curve. In contrast to Fig. 3(b), the overall slope of a trajectory 
now depends on its depth. This important effect, reflected in the different slopes of the 
trajectories in Fig. 4(b), leads to enhanced transport, in the sense that particles at different 
depths have very different histories. Furthermore, the slope (or lateral drift rate) of a 
trajectory is not monotonic in the depth, since it is due to a complicated resonance between 
the lateral motion and downslope topography. 

If the inclination angle 9 is too small, or is such that the furrows are too aligned with the 
incoming flow, then untrapped trajectories disappear but the flow also develops dry spots, 
since the fluid doesn't make it over the crests of the bumps: it is flowing in the 'gutters' of 
the substrate. This takes us outside the scope of the present theory and we do not consider 
such cases. 
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B. Chaotic Trajectories 

In Section VIA we looked at trajectories for the simplest type of substrate, those possess- 
ing a translational symmetry. We now want to investigate the potential for creating flows 
that exhibit chaotic trajectories, leading to chaotic advection. 2-4 Breaking the translational 
symmetry of the substrate is essential for the existence of chaotic trajectories, otherwise the 
trajectories obey a two-dimensional autonomous system, as in Section VIA, which cannot 
exhibit chaos. Given that this translational symmetry is broken, there are three modifica- 
tions to the flow that increase the chances of developing chaotic trajectories: 

1. Break discrete symmetries of the substrate; 

2. Make the fluid deeper (but still thin); 

3. Make the inclination angle 9 smaller. 

The first point removes barriers to transport associated with symmetries. 21-25 The second 
point allows particles more vertical freedom to break the two-dimensional constraint that 
leads to integrable motion, by bringing in first-order terms in the velocity field. The third 
point allows the particles more time to move laterally during one periodic length of the 
domain by decreasing the flow velocity. Here, by laterally we mean a direction perpendicular 
to the mean flow. 

As an illustration we will use the substrate with height function given by (V.3). The first 
term in (V.3) breaks the translational symmetry, and the second term is a perturbation that 
breaks the discrete symmetry x 2 — > x 2 + |. Thus, the substrate with those parameter values 
has a good chance of exhibiting chaos, according to our first point above. 

Consider for an inclination angle 9 = 0.2 for the substrate (V.3), with a corner thick- 
ness e = 0.06. This means that the fluid is reasonably thick, so our second point is satisfied. 
But the Poincare section in Fig. 5(a) clearly shows that no chaos is present, since the phase 
space is foliated by one-dimensional curves. (Though very small regions of chaos could in 
principle be present, they are clearly not important here.) Figure 5(b) shows four typical 
trajectories, taken at about the same initial x 2 but at different depths y. The three un- 
trapped trajectories undergo 'near-trapping,' in the sense that they appear to flirt with the 
idea of entering the trapped regions. This motion is periodic, but the period depends on 
the trajectory. There is thus a depth dependence as for the case in Fig. 4. 

Now we apply our third modification from the list above and decrease the inclination 
angle to 9 = 0.1. This gives the thickness field depicted in Fig. 2. Figure 6(a) shows a 
Poincare section for the resulting flow. The phase space clearly exhibits considerable chaotic 
behavior, as evidenced by the Poincare section not being foliated by one-dimensional surfaces 
[compare to Fig. 5(a)]. There are still a few regions of regular behavior, indicated by the 
four visible islands, which correspond to islands in Fig. 5(a). In Figure 6(b) six trajectories 
are plotted: the bottom three are trapped in the islands. The top three undergo transitions 
between long untrapped motions, where the trajectories go diagonally across the substrate, 
and short trapped motions. This trapped-untrapped behavior is typical of chaotic dynamics 
in the neighborhood of islands. 26 ' 27 
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FIG. 5: For the substrate shape (V.3) with corner-thickness e = 0.06 and inclination angle 6 = 0.1: 
(a) Poincare section at x 1 = 0, showing no apparent chaos; (b) Four typical trajectories, showing 
a trapped trajectory and three trajectories that undergo regular 'near-trapping' events. Note that 
the overall slopes of the trajectories are different, and are not monotonic with depth. 




2 1 
X X 

(a) (b) 

FIG. 6: For same parameters as the flow used for Fig. 5 but with a shallower inclination angle 9 = 
0.1: (a) Poincare section at x 1 = 0, showing considerable chaos; (b) Six typical trajectories, showing 
three trapped trajectories and three trajectories that undergo trapping events and flights typical 
of chaotic dynamics. 
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VII. DISCUSSION 

We have derived an equation describing the gravity-driven motion of a thin viscous flow 
on a curved substrate. This is the same equation as in RRS 8 expressed in nonorthogonal 
coordinates, and in addition including second-order correction terms to satisfy exactly the 
kinematic boundary condition at the free surface. This is essential for particle advection, 
since otherwise the small errors at the surface are enough to cause the particles to escape 
the domain. 

We then investigated the character of fluid trajectories for four cases of increasing com- 
plexity. We observe a hierarchy of complex motions: 

• For a substrate with translational symmetry at a large inclination angle, the flow is 
very thin and particle moves in uniform layers, as in Fig. 3(a). The lateral motion of 
the particles is essentially independent of their depth [Fig. 3(b)]. 

• For the same substrate with a smaller inclination angle, the flow becomes thicker over 
furrows, and zones of recirculation form. These are seen as islands in Fig. 4(a). Two 
new effects are visible in the particles trajectories [Fig. 4(b)]: some are trapped inside 
the rolls, and others exhibit different rates of lateral displacement in x 2 . This means 
that this flow has much better transport properties since particles initially near each 
other will separate rapidly (but not exponentially since there is no chaos). 

• If we break the translational symmetry and use a fully two-dimensional surface such 
as (V.3), we observe the formation of more islands (four visible in total) of trapped 
trajectories [Fig. 5(a)]. The trajectories 'bounce off' the various islands causing com- 
plicated, but regular, oscillations [Fig. 5(b)]. 

• Finally, decreasing the inclination angle of the previous substrate preserves the same 
four islands [Fig. 6(a)], but now the trajectories undergo a succession of random 
trapped and untrapped events, analogous to Levy flights observed in experiments 
in a rotating tank [Fig. 6(b)]. 26,27 This is the most effective transport achieved here, 
since particles initially near each other will separate exponentially, unless they are 
both inside the same island. 

We used three modifications to the flow to increase the complexity of particle trajectories: 
1. Breaking discrete symmetries of the substrate; 2. Making the fluid deeper; 3. Making 
the inclination angle smaller. Of the three, we expect the first two to have a similar effect 
on any other type of thin flow, say one driven by Marangoni stresses or shear at the free 
surface. The last simply says that decreasing the mean velocity of the flow can make lateral 
motions more important, heightening the chance of observing complex dynamics. 

The lateral motion (that is, perpendicular to the mean flow direction) of particles is fairly 
small. We have made little effort to optimize this, preferring to focus our attention on test 
cases as a first object of study. Experiments would be desirable here to see if the observed 
lateral motion is of the right order of magnitude. 

It is tempting to try to compute numerically an effective diffusion coefficient by simulating 
large numbers of particle trajectories. However, the numerical codes used here are too slow 
to achieve the high statistics needed to measure these coefficients with confidence. Rather, 
we hope to use techniques from homogenization theory to calculate the large-scale transport 
properties. 28 
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The generalized coordinate setting introduced in Section II is applicable to any set of 
dynamical equations, and could be used to model flows over complicated surbstrates in a 
fairly straightforward manner. It also allows for a time-dependent substrate, such as when 
it grows by solidification, 9 though we have not treated such effects here. 

Finally, the substrates used here were periodic in both directions, but it would be useful 
to quantify transport for random landscapes, consisting for instance of imperfections on the 
surface of a material. 
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